# This file creates Table SI.3

# install.packages("xtable")
# install.packages("ggplot2")
# install.packages("dplyr")

library(xtable)
library(ggplot2)
library(dplyr)

# read output tables 
outtab_mom <- read.table("dataan/data/mother_firsttwo_res.txt", 
                         header = FALSE, skip = 1)

outtab_dad <- read.table("dataan/data/father_firsttwo_res.txt", 
                         header = FALSE, skip = 1)

main_tab_mom <- 
  matrix(NA, nrow = 5, ncol = 3)

colnames(main_tab_mom) <-
  c("No. of children 2009 Local", "No. of children 2013 Local", "No. of children 2014 European")

rownames(main_tab_mom) <- 
  c("Two first children",
    "same sex",
    "3",
    "No. of children if two first children opposite sex",
    "N")

for(i in 1:3){
  main_tab_mom[1, i] <- round(outtab_mom[2*i + 12, 2], 3)
  main_tab_mom[2, i] <- paste("(",round(outtab_mom[2*i + 12, 3], 3),")", sep ="")
  main_tab_mom[3, i] <- paste("[",
                              round(outtab_mom[2*i + 12, 3]*qnorm(0.025) + outtab_mom[2*i + 12, 2], 3),
                              "; ",
                              round(outtab_mom[2*i + 12, 3]*qnorm(0.975) + outtab_mom[2*i + 12, 2], 3),
                              "]", sep ="")
  main_tab_mom[4, i] <- round(outtab_mom[2*i + 11, 2], 3)
  main_tab_mom[5, i] <- sum(outtab_mom[c(4*i - 3, 4*i - 2), 6])
}



# Set table up for fathers

main_tab_dad <- 
  matrix(NA, nrow = 5, ncol = 3)

colnames(main_tab_dad) <-
  c("No. of children 2009 Local", "No. of children 2013 Local", "No. of children 2014 European")

rownames(main_tab_dad) <- 
  c("Two first children",
    "same sex",
    "3",
    "No. of children if two first children opposite sex",
    "N")

for(i in 1:3){
  main_tab_dad[1, i] <- round(outtab_dad[2*i + 12, 2], 3)
  main_tab_dad[2, i] <- paste("(",round(outtab_dad[2*i + 12, 3], 3),")", sep ="")
  main_tab_dad[3, i] <- paste("[",
                              round(outtab_dad[2*i + 12, 3]*qnorm(0.025) + outtab_dad[2*i + 12, 2], 3),
                              "; ",
                              round(outtab_dad[2*i + 12, 3]*qnorm(0.975) + outtab_dad[2*i + 12, 2], 3),
                              "]", sep ="")
  main_tab_dad[4, i] <- round(outtab_dad[2*i + 11, 2], 3)
  main_tab_dad[5, i] <- sum(outtab_dad[c(4*i - 3, 4*i - 2), 6])
}


## Effects in table
xtable(main_tab_mom)
xtable(main_tab_dad)